Intestinal mTORC2 signaling alters sensory driven traits

Fig 1 Intestinal mTORC2 prevents dauer formation

1B

strains<-c("N2", "rict-1(mg360)", "rict-1(ft7)", "sgk-1(ok538)", "akt-1(mg306)", "akt-2", "pkc-1", "pkc-2")
foods <- "OP50"

TORC2<-read.csv(file.path(pathname, "data", "1B_1C_rict-1_TORC2.csv"), header=TRUE) %>% format_dauer_pd() #including partial/pd as non-dauer due to excess zeros. 

##glmm
#glmm <- TORC2 %>% dauer_pd_glmm() #including partial/pd as non-dauer due to excess zeros. 
#lm
lm <- TORC2 %>% dauer_ANOVA()
#stan
stan.glmm <- TORC2 %>% dauer_stan_glmm()

contrasts<-dauergut::dunnett_contrasts(lm, ref.index = 1, "genotype")

## [1] "interaction term not indicated"

mixed<-stan.glmm %>% dauer_getStan_CIs()
plot.contrasts<-c("",contrasts$prange[1:7])

(p<-dauergut::plot_CIs(TORC2, title='mTORC2 prevents high temperature dauer formation', plot.contrasts, ypos = 1.075, type = "dauer"))

1B. mTORC2 controls high temperature-induced dauer formation. Animals were grown 60hrs post egg lay at 27º, proportion of dauers in each strain is shown. Arrested L3 and/or partial/post dauers are included as non-dauers due to excess arrest in rict-1(ft7). Increased dauer formation occurs in mTORC2-specific rict-1 mutants, as well as in mutants of mTORC2 targets sgk-1 and akt-1. Box and scatter plot show raw data. Bayesian 90% (grey) and 75% (thick-grey) credible intervals are shown on the right (see methods). All p-values reflect ANOVA with Dunnett post-hoc comparision to N2. n≥5 experiments over 3 independent days.

contrast	estimate	SE	df	t.ratio	p.value	prange
rict-1(mg360) - N2	0.6985076	0.0364356	107	19.1710147	0.0000000	p<0.001
rict-1(ft7) - N2	0.6412468	0.0467126	107	13.7274777	0.0000000	p<0.001
sgk-1(ok538) - N2	0.3024057	0.0467126	107	6.4737442	0.0000000	p<0.001
akt-1(mg306) - N2	0.7173955	0.0499946	107	14.3494505	0.0000000	p<0.001
akt-2 - N2	-0.0033547	0.0520678	107	-0.0644299	1.0000000	p~1
pkc-1 - N2	-0.2269687	0.0731846	107	-3.1013152	0.0164674	p~0.016
pkc-2 - N2	-0.2014919	0.0499946	107	-4.0302714	0.0007241	p<0.001

1C

strains<-c("N2","rict-1(mg360)",
           "rict-1(mg360); ex[ges1]",
           "rict-1(mg360); ex[gpa4]",
           "rict-1(ft7)",
           "rict-1(ft7); ex[ges1]",
           "rict-1(ft7); ex[elt2]",
           "rict-1(ft7); ex[ifb2]",
           "sgk-1(ok538)",
           "sgk-1(ok538); ex[ges1]",
           "akt-1(mg306)",
           "akt-1(mg306); ex[ges1]")

gutresc<-read.csv(file.path(pathname, "data", "1C_ges-1_rescue_rict_sgk.csv")) %>% format_dauer_pd() %>% 
  mutate(allele = factor(allele, levels = c("WT","mg360","ft7","ok538","mg306")))

gutresc %<>% mutate(adj.pct = case_when(.$pct == 0 ~ 0.01, .$pct == 1 ~ 0.99, TRUE ~ .$pct))

lm <- gutresc %>% dauer_ANOVA()
stan

## function (file, model_name = "anon_model", model_code = "", fit = NA, 
##     data = list(), pars = NA, chains = 4, iter = 2000, warmup = floor(iter/2), 
##     thin = 1, init = "random", seed = sample.int(.Machine$integer.max, 
##         1), algorithm = c("NUTS", "HMC", "Fixed_param"), control = NULL, 
##     sample_file = NULL, diagnostic_file = NULL, save_dso = TRUE, 
##     verbose = FALSE, include = TRUE, cores = getOption("mc.cores", 
##         1L), open_progress = interactive() && !isatty(stdout()) && 
##         !identical(Sys.getenv("RSTUDIO"), "1"), ..., boost_lib = NULL, 
##     eigen_lib = NULL) 
## {
##     dot_arg_names <- names(list(...))
##     is_arg_recognizable(dot_arg_names, c("chain_id", "init_r", 
##         "test_grad", "append_samples", "refresh", "control", 
##         "enable_random_init", "save_warmup", "obfuscate_model_name"), 
##         pre_msg = "passing unknown arguments: ", call. = FALSE)
##     if (!is.list(data) && !is.environment(data) && !is.character(data)) 
##         stop("'data' must be a list, environment, or character vector")
##     if (is(fit, "stanfit")) 
##         sm <- get_stanmodel(fit)
##     else {
##         attr(model_code, "model_name2") <- deparse(substitute(model_code))
##         if (missing(model_name)) 
##             model_name <- NULL
##         sm <- stan_model(file, model_name = model_name, model_code = model_code, 
##             stanc_ret = NULL, boost_lib = boost_lib, eigen_lib = eigen_lib, 
##             save_dso = save_dso, verbose = verbose)
##     }
##     if (is.null(sample_file)) 
##         sample_file <- NA
##     if (is.null(diagnostic_file)) 
##         diagnostic_file <- NA
##     sampling(sm, data, pars, chains, iter, warmup, thin, seed, 
##         init, check_data = TRUE, sample_file = sample_file, diagnostic_file = diagnostic_file, 
##         verbose = verbose, algorithm = match.arg(algorithm), 
##         control = control, check_unknown_args = FALSE, cores = cores, 
##         open_progress = open_progress, include = include, ...)
## }
## <environment: namespace:rstan>

stan.glmm <- gutresc[!gutresc$pct %in% c(0,1), ] %>% dauer_stan_glmm()

contrasts<-dauergut::tukey_contrasts(lm, "genotype") 
mixed<-stan.glmm %>% dauer_getStan_CIs()
plot.contrasts<-c("",contrasts$prange[1], "", "", contrasts$prange[4], "","","", contrasts$prange[8],"",contrasts$prange[10],"")
plot.contrasts.2<-c("", "",contrasts$prange[12:13],"",contrasts$prange[39:41],"", contrasts$prange[61],"", contrasts$prange[66])

(p<-dauergut::plot_CIs(gutresc, title='mTORC2 components act in the intestine to regulate 
dauer formation', plot.contrasts, plot.contrasts.2=plot.contrasts.2, ypos = 1.20,type = "dauer", offset = 0))

1C. mTORC2 acts in the intestine to inhibit high-temperature dauer formation. Animals were grown 60hrs post egg lay at 27º, proportion of dauers in each strain is shown. Increased dauer formation of rict-1 and sgk-1 mutants is rescued by intestinal-specific expression. akt-1 mutants were not rescued using this promoter. Box and scatter plot show raw data. Bayesian 90% (grey) and 50% (red) credible intervals are shown on the right (see methods). All p-values reflect ANOVA with Dunnett post-hoc comparision to N2. p-values in red indicate rescue comparisons to respective control mutants. n≥6 experiments over 3 independent days.

contrast	estimate	SE	df	t.ratio	p.value	prange
N2 - rict-1(mg360)	-0.6909880	0.0518355	100	-13.3303901	0.0000000	p<0.001
N2 - rict-1(mg360); ex[ges1]	0.1241772	0.0476694	100	2.6049685	0.2828932	p~0.28
N2 - rict-1(mg360); ex[gpa4]	-0.6605064	0.0460611	100	-14.3397941	0.0000000	p<0.001
N2 - rict-1(ft7)	-0.6599220	0.0476694	100	-13.8437300	0.0000000	p<0.001
N2 - rict-1(ft7); ex[ges1]	0.0305997	0.0518355	100	0.5903222	0.9999818	p~1
N2 - rict-1(ft7); ex[elt2]	-0.1258822	0.0518355	100	-2.4284924	0.3871594	p~0.39
N2 - rict-1(ft7); ex[ifb2]	-0.0552211	0.0581192	100	-0.9501352	0.9982031	p~1
N2 - sgk-1(ok538)	-0.4428609	0.0581192	100	-7.6198664	0.0000000	p<0.001
N2 - sgk-1(ok538); ex[ges1]	0.0831857	0.0581192	100	1.4312943	0.9512955	p~0.95
N2 - akt-1(mg306)	-0.7145009	0.0581192	100	-12.2937048	0.0000000	p<0.001
N2 - akt-1(mg306); ex[ges1]	-0.7176764	0.0581192	100	-12.3483425	0.0000000	p<0.001
rict-1(mg360) - rict-1(mg360); ex[ges1]	0.8151652	0.0610818	100	13.3454703	0.0000000	p<0.001
rict-1(mg360) - rict-1(mg360); ex[gpa4]	0.0304817	0.0598351	100	0.5094278	0.9999960	p~1
rict-1(mg360) - rict-1(ft7)	0.0310660	0.0610818	100	0.5085968	0.9999961	p~1
rict-1(mg360) - rict-1(ft7); ex[ges1]	0.7215877	0.0643859	100	11.2072381	0.0000000	p<0.001
rict-1(mg360) - rict-1(ft7); ex[elt2]	0.5651058	0.0643859	100	8.7768615	0.0000000	p<0.001
rict-1(mg360) - rict-1(ft7); ex[ifb2]	0.6357669	0.0695447	100	9.1418477	0.0000000	p<0.001
rict-1(mg360) - sgk-1(ok538)	0.2481271	0.0695447	100	3.5678807	0.0250011	p~0.025
rict-1(mg360) - sgk-1(ok538); ex[ges1]	0.7741738	0.0695447	100	11.1320344	0.0000000	p<0.001
rict-1(mg360) - akt-1(mg306)	-0.0235128	0.0695447	100	-0.3380969	0.9999999	p~1
rict-1(mg360) - akt-1(mg306); ex[ges1]	-0.0266883	0.0695447	100	-0.3837582	0.9999998	p~1
rict-1(mg360); ex[ges1] - rict-1(mg360); ex[gpa4]	-0.7846836	0.0562644	100	-13.9463561	0.0000000	p<0.001
rict-1(mg360); ex[ges1] - rict-1(ft7)	-0.7840992	0.0575885	100	-13.6155604	0.0000000	p<0.001
rict-1(mg360); ex[ges1] - rict-1(ft7); ex[ges1]	-0.0935776	0.0610818	100	-1.5320041	0.9238594	p~0.92
rict-1(mg360); ex[ges1] - rict-1(ft7); ex[elt2]	-0.2500595	0.0610818	100	-4.0938460	0.0044542	p~0.004
rict-1(mg360); ex[ges1] - rict-1(ft7); ex[ifb2]	-0.1793984	0.0664974	100	-2.6978240	0.2353514	p~0.24
rict-1(mg360); ex[ges1] - sgk-1(ok538)	-0.5670381	0.0664974	100	-8.5272183	0.0000000	p<0.001
rict-1(mg360); ex[ges1] - sgk-1(ok538); ex[ges1]	-0.0409915	0.0664974	100	-0.6164370	0.9999715	p~1
rict-1(mg360); ex[ges1] - akt-1(mg306)	-0.8386781	0.0664974	100	-12.6121875	0.0000000	p<0.001
rict-1(mg360); ex[ges1] - akt-1(mg306); ex[ges1]	-0.8418536	0.0664974	100	-12.6599413	0.0000000	p<0.001
rict-1(mg360); ex[gpa4] - rict-1(ft7)	0.0005843	0.0562644	100	0.0103857	1.0000000	p~1
rict-1(mg360); ex[gpa4] - rict-1(ft7); ex[ges1]	0.6911060	0.0598351	100	11.5501788	0.0000000	p<0.001
rict-1(mg360); ex[gpa4] - rict-1(ft7); ex[elt2]	0.5346241	0.0598351	100	8.9349594	0.0000000	p<0.001
rict-1(mg360); ex[gpa4] - rict-1(ft7); ex[ifb2]	0.6052852	0.0653541	100	9.2616222	0.0000000	p<0.001
rict-1(mg360); ex[gpa4] - sgk-1(ok538)	0.2176455	0.0653541	100	3.3302483	0.0504653	p~0.05
rict-1(mg360); ex[gpa4] - sgk-1(ok538); ex[ges1]	0.7436921	0.0653541	100	11.3794210	0.0000000	p<0.001
rict-1(mg360); ex[gpa4] - akt-1(mg306)	-0.0539945	0.0653541	100	-0.8261836	0.9995050	p~1
rict-1(mg360); ex[gpa4] - akt-1(mg306); ex[ges1]	-0.0571700	0.0653541	100	-0.8747727	0.9991556	p~1
rict-1(ft7) - rict-1(ft7); ex[ges1]	0.6905217	0.0610818	100	11.3048694	0.0000000	p<0.001
rict-1(ft7) - rict-1(ft7); ex[elt2]	0.5340398	0.0610818	100	8.7430275	0.0000000	p<0.001
rict-1(ft7) - rict-1(ft7); ex[ifb2]	0.6047009	0.0664974	100	9.0935972	0.0000000	p<0.001
rict-1(ft7) - sgk-1(ok538)	0.2170611	0.0664974	100	3.2642030	0.0606100	p~0.061
rict-1(ft7) - sgk-1(ok538); ex[ges1]	0.7431078	0.0664974	100	11.1749842	0.0000000	p<0.001
rict-1(ft7) - akt-1(mg306)	-0.0545788	0.0664974	100	-0.8207663	0.9995353	p~1
rict-1(ft7) - akt-1(mg306); ex[ges1]	-0.0577544	0.0664974	100	-0.8685200	0.9992094	p~1
rict-1(ft7); ex[ges1] - rict-1(ft7); ex[elt2]	-0.1564819	0.0643859	100	-2.4303767	0.3846168	p~0.38
rict-1(ft7); ex[ges1] - rict-1(ft7); ex[ifb2]	-0.0858208	0.0695447	100	-1.2340386	0.9837202	p~0.98
rict-1(ft7); ex[ges1] - sgk-1(ok538)	-0.4734606	0.0695447	100	-6.8080056	0.0000000	p<0.001
rict-1(ft7); ex[ges1] - sgk-1(ok538); ex[ges1]	0.0525861	0.0695447	100	0.7561481	0.9997887	p~1
rict-1(ft7); ex[ges1] - akt-1(mg306)	-0.7451005	0.0695447	100	-10.7139833	0.0000000	p<0.001
rict-1(ft7); ex[ges1] - akt-1(mg306); ex[ges1]	-0.7482760	0.0695447	100	-10.7596446	0.0000000	p<0.001
rict-1(ft7); ex[elt2] - rict-1(ft7); ex[ifb2]	0.0706611	0.0695447	100	1.0160530	0.9967308	p~1
rict-1(ft7); ex[elt2] - sgk-1(ok538)	-0.3169787	0.0695447	100	-4.5579140	0.0007995	p<0.001
rict-1(ft7); ex[elt2] - sgk-1(ok538); ex[ges1]	0.2090680	0.0695447	100	3.0062396	0.1172395	p~0.117
rict-1(ft7); ex[elt2] - akt-1(mg306)	-0.5886186	0.0695447	100	-8.4638917	0.0000000	p<0.001
rict-1(ft7); ex[elt2] - akt-1(mg306); ex[ges1]	-0.5917941	0.0695447	100	-8.5095530	0.0000000	p<0.001
rict-1(ft7); ex[ifb2] - sgk-1(ok538)	-0.3876397	0.0743464	100	-5.2139687	0.0000517	p<0.001
rict-1(ft7); ex[ifb2] - sgk-1(ok538); ex[ges1]	0.1384069	0.0743464	100	1.8616491	0.7706380	p~0.77
rict-1(ft7); ex[ifb2] - akt-1(mg306)	-0.6592797	0.0743464	100	-8.8676763	0.0000000	p<0.001
rict-1(ft7); ex[ifb2] - akt-1(mg306); ex[ges1]	-0.6624552	0.0743464	100	-8.9103885	0.0000000	p<0.001
sgk-1(ok538) - sgk-1(ok538); ex[ges1]	0.5260466	0.0743464	100	7.0756178	0.0000000	p<0.001
sgk-1(ok538) - akt-1(mg306)	-0.2716400	0.0743464	100	-3.6537076	0.0192617	p~0.019
sgk-1(ok538) - akt-1(mg306); ex[ges1]	-0.2748155	0.0743464	100	-3.6964198	0.0168881	p~0.017
sgk-1(ok538); ex[ges1] - akt-1(mg306)	-0.7976866	0.0743464	100	-10.7293254	0.0000000	p<0.001
sgk-1(ok538); ex[ges1] - akt-1(mg306); ex[ges1]	-0.8008621	0.0743464	100	-10.7720376	0.0000000	p<0.001
akt-1(mg306) - akt-1(mg306); ex[ges1]	-0.0031755	0.0743464	100	-0.0427122	1.0000000	p~1

Fig 2 Intestinal mTORC2 regulates sensory neuron state

2A_p1

Image of daf-7 expression

2A

strains<-c("N2","mg360","ft7","ft7_ges1resc")
dates<-c("10_22_16", "11_9_16", "11_23_16") #dropped "12_6_16" due to missingness
d7GFP<-read.csv("data/2A_daf7GFP.csv") %>% filter(mean!=4095 & genotype %in% strains & date %in% dates & temp == "27" & food == "OP50") %>% mutate(genotype = factor(genotype, levels=strains), ID = as.character(ID)) %>%
  separate(ID, c("ID.A", "ID.B"), sep = ":", extra = "drop") %>% 
  mutate(genoID = paste(date, genotype, ID.A, sep = ":"), cell.norm = mean) #genoID is animal, 2 cells per animal measured.

df <- d7GFP %>% group_by(date, genotype, genoID) %>% summarise(cell.norm = mean(cell.norm)) # take mean of each worm

lmm1<-lmer(cell.norm ~ genotype + (1|date), data = d7GFP)
library(DHARMa)
#plotSimulatedResiduals(simulationOutput = simulateResiduals(lmm1, n=1000))
# plot of residuals showed heteroskedastic data (increasing by mean), log transforming the response corrected this.
library(lsmeans)
# log.tran<-make.tran(type="genlog", param = c(0.5,10))
log.tran<-make.tran(type="genlog", param = c(0,10))
lmm2<-with(log.tran, lmer(linkfun(cell.norm) ~ genotype + (1|date) + (1|genotype:date), data = df))
lm <- with(log.tran, lm(linkfun(cell.norm) ~ genotype, data = df))
stanlmer <- with(log.tran, stan_lmer(linkfun(cell.norm) ~ genotype + (1|date) + (1:genotype:date), data = df))

#compare ft7 to rescue H0, ft7 = ft7_resc
rescue.test <- d7GFP %>% subset(genotype %in% c("ft7", "ft7_ges1resc")) %$% with(log.tran, t.test(linkfun(cell.norm) ~ genotype))
rescue.p <- data.frame(p.value = rescue.test$p.value) %>% dauergut::prange()

#H0, all genotypes = N2
contrasts<-dauergut::dunnett_contrasts(lm, ref.index = 1, "genotype")

## [1] "interaction term not indicated"

mixed <- stanlmer %>% getStan_CIs_log(base = 10)

plot.contrasts<-c("",contrasts$prange[1:2],"")
plot.contrasts.2 <- c("", "", "", rescue.p$prange)

dauergut::plot_CIs(df=df, title='daf-7 expression in ASI at 27º is
controlled by intestinal rict-1', plot.contrasts=plot.contrasts, plot.contrasts.2 = plot.contrasts.2, ypos = 4700, offset = 0, type = "GFP")

2A. Intesinal rict-1 regulates daf7GFP levels in ASI neurons. L1 animals were grown 18hrs post egg lay at 27º. Daf-7::GFP fluorescence quantification shows decreased expression in rict-1 mutants. Intestinal expression of rict-1 using the intestinal ges-1 promoter rescues daf-7::GFP expression. Box and scatter plot show raw data. Bayesian 90% (grey) and 50% (red) credible intervals are shown on the right (see methods). For rict-1(mg360) and rict-1(ft7), p<0.001 and p~0.003 compared to N2, respectively. For ges-1 rescue, p<0.001 compared to rict-1(ft7). n=3 independent experiments.

knitr::kable(contrasts, caption="Pairwise comparisons from mixed-effects model (Tukey)")

Pairwise comparisons from mixed-effects model (Tukey)
contrast	estimate	SE	df	t.ratio	p.value	prange
mg360 - N2	-0.3044202	0.0444226	122	-6.8528183	0.0000000	p<0.001
ft7 - N2	-0.1504361	0.0447683	122	-3.3603231	0.0030178	p~0.003
ft7_ges1resc - N2	0.0151539	0.0463634	122	0.3268511	0.9772211	p~0.98

knitr::kable(mixed[,c(6,1:5)], caption = "Bayesian credible intervals")

Bayesian credible intervals
genotype	mean	lower.CL	upper.CL	lower.25	upper.75
N2	2328.037	1808.907	3758.711	2029.522	2638.144
mg360	1164.326	956.080	1412.707	1038.259	1308.172
ft7	1648.532	1360.950	2003.864	1469.812	1849.672
ft7_ges1resc	2403.224	1946.362	2970.160	2120.776	2720.512

2B

daf7FISH<-read.csv("data/2B_daf7FISH_4frames.csv") %>%
  separate(Label, c("method","group","sample"),sep = c("-"), extra = "drop") %>%
  separate(sample, c("sample", "type"),sep =":", extra = "drop") %>%
  mutate(genotype = data.frame(do.call(rbind, strsplit(as.vector(group), split = "daf7_mRNA_L1_27d_")))[,2]) %>%
  mutate(ID = interaction(genotype, sample)) %>%
  mutate(genotype = factor(genotype, levels = c("N2", "ft7", "mg360", "mg360resc")))

#measured background fluorescence of each of 4 frames to normalize the sum projection of the neuron
#merge to add mean background value to each row:
background <- filter(daf7FISH, type == "background") %>% dplyr::select(ID, mean)
colnames(background) <- c("ID", "background")

cells<-dplyr::filter(daf7FISH, type == "cell") %>% dplyr::select(ID, sample, mean, genotype)
colnames(cells) <- c("ID", "sample", "cell.mean", "genotype")

daf7FISH <- merge(cells, background,by="ID") %>%
  mutate(cell.norm = cell.mean - mean(background), 
         cell.diffnorm = cell.mean - background,
         cell.ratioNorm = cell.mean / background)

#simple lm
lm<-lm(data=daf7FISH, formula=cell.norm~genotype)
#check for outliers
daf7FISH <- dauergut::flag_outliers(lm, daf7FISH, threshold = 4, noplot = TRUE)

lm.0<-lm(data=daf7FISH[daf7FISH$outlier.status == FALSE,], formula=cell.norm~1)
lm.1<-lm(data=daf7FISH[daf7FISH$outlier.status == FALSE,], formula=cell.norm~genotype)
#anova(lm.0, lm.1) p ~ 0.013

stanlm <- stan_glm(cell.norm ~ genotype, data = daf7FISH[daf7FISH$outlier.status == FALSE,])

strains <- levels(daf7FISH$genotype)
contrasts <- dauergut::dunnett_contrasts(lm.1,ref.index = 1,"genotype")

## [1] "interaction term not indicated"

mixed <- stanlm %>% getStan_CIs()
rescue.test <- daf7FISH[daf7FISH$outlier.status == FALSE,] %>% subset(genotype %in% c("mg360", "mg360resc")) %$% t.test(cell.norm ~ genotype)
rescue.p <- data.frame(p.value = rescue.test$p.value) %>% dauergut::prange()
#mixed <-dauergut::MM_CI_trans_dunnett(lm, strains)


plot.contrasts <- c("", contrasts$prange[1:2],"")
plot.contrasts.2 <- c("", "", "", rescue.p$prange)

dauergut::plot_CIs(daf7FISH[daf7FISH$outlier.status == FALSE,], title = "daf7 mRNA is reduced in rict-1 mutants", plot.contrasts = plot.contrasts, plot.contrasts.2 = plot.contrasts.2, ypos = 150, offset = 0, type = "expression")

2C. Daf-7 mRNA is reduced in rict-1 mutants. Animals were grown 18hrs post egg lay at 27º, then fixed prior to FISH. Intestinal expression of rict-1 using the intestinal ges-1 promoter rescues daf-7 expression. Box and scatter plot show raw data. Bayesian 90% (grey) and 50% (red) credible intervals are shown on the right (see methods). For rict-1(mg360) and rict-1(ft7), p~0.04 and p~0.032 compared to N2, respectively. For ges-1 rescue, p~0.038 compared to rict-1(ft7). n=2 independent experiments, pooled data. 2 outlier data points were removed in total.

2C

rescue of rict-1 with daf-7 and daf-28

strains<-c("N2","rict-1(ft7)","rict-1(ft7); ex[ASI::daf7]","rict-1(ft7); ex[ASI::daf28]","rict-1(ft7); ex[ASJ::daf28]")
daf7supp<-read.csv(file.path(pathname, "data/2C_daf7_daf28_suppression.csv")) %>% format_dauer_pd()

# daf7supp$genotype<- factor(daf7supp$genotype, levels = strains)
# daf7supp$pct<-as.numeric(paste(daf7supp$dauer/(daf7supp$dauer+daf7supp$pd + daf7supp$non))) #omitted arrest
# daf7supp$non.dauer<-as.numeric(paste(daf7supp$pd+daf7supp$non))

lm <- daf7supp %>% dauer_ANOVA()
#stan
stan.glmm <- daf7supp %>% dauer_stan_glmm()

contrasts<-dauergut::tukey_contrasts(lm, "genotype")
mixed<-stan.glmm %>% dauer_getStan_CIs()
plot.contrasts<-c("",contrasts$prange[1],"","","")
plot.contrasts.2<-c("", "", contrasts$prange[5:7]) #for rescue vs rict-1

(p<-dauergut::plot_CIs(daf7supp, title='daf-7 expression in ASI suppresses rict-1 dauer phenotype', plot.contrasts, plot.contrasts.2, ypos = 1.075, offset = 0, type = "dauer"))

2D

#days = as.factor(8:10) # days containing N2 and ft7 27º data
days = as.factor(8:10)
strains = c("N2", "ft7")
foods = "OP50"
d28<-read.csv('data/2D_3B_daf-28_GFP_rict-1.csv') %>%
separate(ID, c("ID.A", "ID.B"), sep = ":", extra = "drop") %>%
  subset(mean!=4095 & food == foods & temp == "27" & genotype %in% strains & pheromone == 0) %>% 
  mutate(day = as.factor(day), genotype = factor(genotype, levels = strains),
         genoID = interaction(date,genotype, ID.A))

df <- d28 %>% subset(day %in% days) %>% group_by(date,genotype,neuron,genoID) %>% summarise(cell.norm = mean(mean)) %>%
  data.frame() %>% mutate(group.id = interaction(genotype, neuron)) # take mean of each worm

#fit simple linmod to est outliers:
lm1 <- with(log.tran, lm(linkfun(cell.norm) ~ neuron * genotype, data = df))
df %<>% flag_outliers(lin.mod=lm1, threshold = 4, df = .)

#df.ASI <- df %>% dplyr::filter(neuron == "ASI")
#df.ASJ <- df %>% dplyr::filter(neuron == "ASJ")
# linear MM
library(lsmeans)
log.tran<-make.tran(type="genlog", param = c(0,10)) #make log10 transformation
#lmm.ASI <- with(log.tran, lmer(linkfun(cell.norm) ~ neuron * genotype + (1|neuron:date) + (1|genotype:neuron:date), data = df.ASI[df$outlier.status==FALSE,]))
#lmm.ASJ <- with(log.tran, lmer(linkfun(cell.norm) ~ neuron * genotype + (1|neuron:date) + (1|genotype:neuron:date), data = df.ASJ[df$outlier.status==FALSE,]))
lm <- with(log.tran, lm(linkfun(cell.norm) ~ neuron * genotype, data = df[df$outlier.status==FALSE,]))
stanlmer <- with(log.tran, stan_lmer(linkfun(cell.norm) ~ 1 + group.id + (1|date) + (1:group.id:date), data = df))


#plotSimulatedResiduals(simulationOutput = simulateResiduals(lm, n=1000)) log + outlier removal fixes qq-plot



#get confidence intervals:
# fixed.names <- attributes(summary(brms.lmm)$fixed)$dimnames[[1]] # get names of fixed effects in model:
# h0 <- paste(fixed.names[1], "= 0") #N2, ASI
# h1 <- paste(fixed.names[1],"+",fixed.names[2], "= 0") #ft7, ASI
# h2 <- paste(fixed.names[1],"+",fixed.names[3], "= 0") #N2, ASJ
# h3 <- paste(fixed.names[1],"+",fixed.names[2], "+",fixed.names[3], "+", fixed.names[4], "= 0") #ft7, ASJ
# hyp.list <- c(h0, h1, h2, h3)

# list.90<-lapply(hyp.list, function(x) {
#   brms::hypothesis(brms.lmm, x, alpha = 0.1)[1][[1]]
# }) %>% bind_rows() %>% dplyr::select(1:4) %>% 10^. %>% data.frame() %>% mutate(fixed.effect = fixed.names)
# 
# list.50<-lapply(hyp.list, function(x) {
#   brms::hypothesis(brms.lmm, x, alpha = 0.5)[1][[1]]
# }) %>% bind_rows() %>% dplyr::select(1:4) %>% 10^. %>% data.frame() %>% mutate(fixed.effect = fixed.names)

# CIs <- full_join(list.90, list.50) %>% dplyr::select(1,3:4,5:7) %>% 
#   `colnames<-`(c("mean","lower.CL", "upper.CL", "fixed.effect", "lower.25", "upper.75")) %>% 
#   mutate(genotype = rep(strains, 2),
#          neuron = rep(c("ASI","ASJ"), each = 2),
#          x.pos = rep(c(1.3, 2.3), 2))


# #hypothesis tests:
# h0 <- paste(fixed.names[1], ">" , fixed.names[1], " + ", fixed.names[2]) #ft7 less than N2 (ASI)
# h1 <- paste(fixed.names[1],"+",fixed.names[3], ">" ,
#             fixed.names[1],"+",fixed.names[2], "+",fixed.names[3], "+", fixed.names[4]) #ft7 less than N2 (ASJ)
# hyp.list <- c(h0, h1)
# 
# hyps <- lapply(hyp.list, function(x) {
#   name = x
#   brms::hypothesis(brms.lmm, x, alpha = 0.03)
# })




# # 
# stanglmm <- stan_lmer(cell.norm ~ group.id, data = df,
#   link = "log", prior = normal(), prior_intercept = normal(),
#   prior_aux = cauchy(0, 5))

with(df[df$neuron == "ASI" & df$outlier.status == FALSE,], with(log.tran, t.test(cell.norm ~genotype)))

## 
##  Welch Two Sample t-test
## 
## data:  cell.norm by genotype
## t = 2.1727, df = 38.435, p-value = 0.03603
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   30.1321 848.2671
## sample estimates:
##  mean in group N2 mean in group ft7 
##          1664.254          1225.055

with(df[df$neuron == "ASJ" & df$outlier.status == FALSE,], with(log.tran, t.test(cell.norm ~genotype)))

## 
##  Welch Two Sample t-test
## 
## data:  cell.norm by genotype
## t = 5.0497, df = 53.001, p-value = 0.000005574
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  183.8804 426.2096
## sample estimates:
##  mean in group N2 mean in group ft7 
##         1073.3346          768.2896

#mixed <- CIs



contrasts<-summary(lsmeans(lm, pairwise ~ genotype | neuron)$contrasts) %>% dauergut::prange()
mixed <- getStan_CIs_log_neuron(stanlmer,group="genotype",base=10)
plot.contrasts<-c("",contrasts$prange[1],"", contrasts$prange[2])

#### plot ######
p<-(df %>% ggplot(aes(x=genotype, y = cell.norm)) +
  #### plot raw data and quartiles  
          geom_quasirandom(aes(y=cell.norm, pch = date),cex=1,colour = "#339900",
                           width = 0.075,size=0.3,
                           method = 'smiley') +
    guides(pch=FALSE) +
    list(add.median(width = 0.25), add.quartiles()) +
labs(title = "daf-28 is also reduced in rict-1 mutants",
     subtitle = "GFP plot",
           y = "mean GFP intensity", x="genotype") +
  #### plot Bayesian credible intervals
  stat_summary(aes(x=genotype, group = genotype, y = 3500), geom="text", fun.data=box_annot, label=plot.contrasts, size = 6) +
  geom_errorbar(data=mixed, aes(x=x.pos,y=mean, ymin=lower.CL, ymax=upper.CL),
                  width=0,colour ="grey", lwd=0.2) + #90% cred int for 95% one-sided H0
  geom_errorbar(data=mixed, aes(x=x.pos,y=mean, ymin = lower.25, ymax = upper.75),
                  width=0,colour = "darkgrey", lwd = 0.9) + #75% cred interval
      geom_segment(data = mixed, aes(x = x.pos-(0.007*nrow(mixed)), y = mean, xend = x.pos+(0.007*nrow(mixed)), yend = mean), colour = "darkgrey") + 
  stat_summary(aes(x=as.numeric(as.factor(genotype)) + 0.3,y=0,group = genotype),
               fun.data = fun_length, geom = "text",size = 4) +
  facet_grid(.~neuron, switch="both") + theme_my +
  theme(axis.text.y = element_text(size=12),
        axis.text.x = element_text(size = 16),
        axis.title = element_text(size=16),
        panel.spacing.x=unit(2,"lines")))

# knitr::kable(summary(lsmeans(lm3, pairwise ~ genotype | neuron, type = "response"))$contrasts)

S1A

strains<-c("N2","rict-1(mg360)","rict-1;daf-16")
foods <- "OP50"
daf16supp<-read.csv(file.path(pathname, "data/S2_daf16_suppression.csv")) %>% format_dauer_pd()

# daf28supp$genotype<- factor(daf28supp$genotype, levels = strains)
# daf28supp$pct<-as.numeric(paste(daf28supp$dauer/(daf28supp$dauer+daf28supp$pd + daf28supp$non))) #omitted arrest
# daf28supp$non.dauer<-as.numeric(paste(daf28supp$pd+daf28supp$non))
lm <- daf16supp %>% dauer_ANOVA()
#stan
stan.glmm <- daf16supp %>% dauer_stan_glmm()

contrasts<-dauergut::tukey_contrasts(lm, "genotype")
mixed<-stan.glmm %>% dauer_getStan_CIs()
plot.contrasts<-c("",contrasts$prange[1],"")
plot.contrasts.2<-c("", "", contrasts$prange[3]) #for rescue vs rict-1

(p<-dauergut::plot_CIs(daf16supp, title='rict-1 acts in the insulin pathway to inhibit dauer formation', plot.contrasts, plot.contrasts.2, ypos = 1.075, offset = 0, type = "dauer"))

S1BC

strains <- c("N2", "rict-1(ft7)")
foods <- "OP50"
rictC3 <- read.csv(file.path(pathname, 'data/S3C_rict1Pheromone.csv')) %>% format_dauer() %>% mutate(plateID = interaction(plateID,concentration.uM.))

glm <- glm(data = rictC3, cbind(dauer, non) ~ concentration.uM. * genotype, family = binomial)
glmm <- glmer(data = rictC3, cbind(dauer, non) ~ concentration.uM. * genotype + (1|plateID), family = binomial)
newdata = data.frame(genotype = rep(strains, each = 241), concentration.uM. = seq(0,2.4, by = 0.01))
predictions <- predict(glm, newdata = newdata, type = "response", se.fit = TRUE)
#glmm.predict <- predict(glmm,newdata=newdata,type = "response", re.form = NA) #re.form = ~0) 
newdata1 <- cbind(newdata, predictions)
newdata1 %<>% mutate(lower = (fit - se.fit), upper = (fit + se.fit), pct = fit)
#newdata2 <- cbind(newdata, glmm.predict) %>% mutate(pct = glmm.predict)
#newdata3 <- effects::effect(newdata, glmm)


(p<-rictC3 %>% ggplot(aes(x=concentration.uM.,y = pct)) +
  geom_ribbon(data = newdata1,aes(ymin=lower, ymax=upper,fill=genotype), alpha=0.3) + 
  geom_line(data = newdata1,aes(x = concentration.uM., y = pct, colour = genotype)) + 
  geom_point(aes(y=pct), size = 0.6, alpha = 0.75, width = 0.05) +
  add.median.dauer() +
  #geom_line(data = newdata2,aes(x = concentration.uM., y = pct), colour = "green") + 
  labs(y = "proportion dauer") +
  facet_grid(.~genotype) + #, switch="both") +
  scale_fill_manual(values = c("grey", "lightblue")) +
  scale_colour_manual(values = c("grey", "lightblue")) +
  scale_y_continuous(breaks=c(0,0.25,0.5,0.75, 1.0)) +
  #scale_x_discrete(labels=function(x) sub(" ","\n",x,fixed=TRUE)) +
  #geom_text(data = plot.contrasts.H0, aes(x=2, label = prange, y = 1.075, group = NULL), size = 3) +
  theme_my_classic + 
  theme(#axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_text(size=16),
        axis.title.y = element_text(size =20),
        strip.text.x = element_blank(),
        panel.spacing = unit(2,"lines"))
)

Fig 3 Daf-28, but not daf-7, is strongly reduced by temperature

3A

daf-7 temp curve

strains = "N2"
#strains<-c("N2","ft7")
days<-c("11_23_16", "12_6_16")
#days <- "11_23_16" # 12_6 was likely L2s based on the timing

d7GFP<-read.csv(file.path(pathname, "data/3A_daf7GFP_nopeptone_plates.csv")) %>% filter(mean!=4095 & genotype %in% strains & date %in% days & pheromone == 0 & food == "OP50") %>% mutate(genotype = factor(genotype, levels=strains), ID = as.character(ID)) %>%
  separate(ID, c("ID.A", "ID.B"), sep = ":", extra = "drop") %>% 
  mutate(genoID = as.factor(paste(date, genotype, ID.A, sep = ":")), cell.norm = mean)

df <- d7GFP %>% group_by(genotype, date, neuron, temp, genoID) %>% summarise(cell.norm = mean(cell.norm))

#### plot ####
p<-df %>% ggplot(aes(x=temp, y=cell.norm)) +
  list(add.scatter(),add.median(width = .5),add.quartiles()) + 
  labs(
  title = "daf-7 expression does not decrease at high temperatures",
  x = "temperature (ºC)",
  y = "mean intensity"
) + stat_smooth(method="lm") +
  add.n() +
  theme_classic() + 
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 15),
        strip.text.x = element_text(size=20))
p

lm<- lm(cell.norm ~ temp, data = d7GFP)

3B

tempdays<-c("1", "4", "6", "8", "10")
foods <- "OP50"
strains <- "N2"
d28<-read.csv('data/2D_3B_daf-28_GFP_rict-1.csv') %>%
separate(ID, c("ID.A", "ID.B"), sep = ":", extra = "drop") %>%
  subset(mean!=4095 & food == foods & genotype == strains & pheromone == 0) %>% 
  mutate(day = as.factor(day), genotype = factor(genotype, levels = strains),
         genoID = interaction(date,neuron,temp,genotype,ID.A))

df <- d28 %>% subset(day %in% tempdays) %>% group_by(date,temp,genotype,neuron,genoID) %>% summarise(cell.norm = mean(mean)) %>%
  data.frame() %>% mutate(group.id = interaction(genotype, neuron)) # take mean of each worm

#### plot ####
(p1<-df %>% ggplot(aes(x=temp, y=cell.norm)) +
   list(add.scatter(),add.median(width = .5),add.quartiles()) + 
    theme_classic() +
  facet_grid(.~neuron) +
  scale_color_viridis(discrete = TRUE) +
  labs(title = "daf-28 expression is controlled by temperature",
    x = "temperature",
    y = "mean fluorescence",
    colour = "Temperature") + add.n() + 
    theme(axis.title.x = element_text(size = 20),
          axis.title.y = element_text(size= 20),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 15),
        strip.text.x = element_text(size=20),
        panel.spacing.x=unit(2, "lines")))

#### Linear MM ###3
library(lsmeans)
log.tran<-make.tran(type="genlog", param = c(0.5,10))
lm4<-with(log.tran, lmer(linkfun(mean) ~ neuron * temp + (1|genoID), data = subset(d28, day %in% tempdays & food == "OP50" & mean !=4095 & genotype == "N2")))
# lm5<-with(log.tran, lmer(linkfun(mean) ~ exp(temp) + neuron + (1|genoID), data = subset(d28, day %in% tempdays & food == "OP50" & mean !=4095 & genotype == "N2")))
# 
# 
# contrasts<-summary(lsmeans(lm4, ~ neuron | temp, type = "response"))$contrasts %>% prange()
# mixed<-summary(lsmeans(lm4, pairwise ~ genotype | neuron, type = "response"))$lsmeans
# mixed$x.pos<-(as.numeric(as.factor(mixed$genotype)) + 0.3)
# 
# 
# plot.contrasts<-c("",contrasts$prange[1:2],"", contrasts$prange[4:5])




p1

#scale_colour_manual(values = c("grey","#FDE725FF", "#440154FF")) +

3C

days = as.factor(8:10) # days containing N2 and ft7 27º data
temps = c("25","27")
strains = c("N2", "ft7")
foods = "OP50"
d28<-read.csv('data/2D_3B_daf-28_GFP_rict-1.csv') %>%
separate(ID, c("ID.A", "ID.B"), sep = ":", extra = "drop") %>%
  subset(mean!=4095 & food == foods &
           temp %in% temps &
           genotype %in% strains &
           pheromone == 0 &
           day %in% days) %>%
  mutate(day = as.factor(day), genotype = factor(genotype, levels = strains),
         genoID = interaction(date,genotype,ID.A), temp = factor(temp, levels = temps))

df <- d28 %>% group_by(date,temp,genotype,neuron,genoID) %>% summarise(cell.norm = mean(mean)) %>%
  data.frame() %>% mutate(group.id = interaction(genotype, temp)) # take mean of each worm

#### add normalized to 25 degree data
# means <- df %>% group_by(genotype, neuron, temp) %>% summarize(mean = mean(cell.norm))
# df.WT.ASI <- df %>% dplyr::filter(genotype == "N2" & neuron == "ASI") %>% mutate(norm.mean = cell.norm/as.numeric(means[1,4]))
# df.ft7.ASI <- df %>% dplyr::filter(genotype == "ft7" & neuron == "ASI") %>% mutate(norm.mean = cell.norm/as.numeric(means[5,4]))
# df.WT.ASJ <- df %>% dplyr::filter(genotype == "N2" & neuron == "ASJ") %>% mutate(norm.mean = cell.norm/as.numeric(means[3,4]))
# df.ft7.ASJ <- df %>% dplyr::filter(genotype == "ft7" & neuron == "ASJ") %>% mutate(norm.mean = cell.norm/as.numeric(means[7,4]))
# 
# df <- rbind(df.WT.ASI,df.ft7.ASI,df.WT.ASJ,df.ft7.ASJ)

##### model specification ####
# lm1 <- with(log.tran, lm(linkfun(cell.norm) ~ neuron, data = df))
# lm2 <- with(log.tran, lm(linkfun(cell.norm) ~ neuron + temp, data = df))
# lm3 <- with(log.tran, lm(linkfun(cell.norm) ~ neuron * temp, data = df))
# lm4 <- with(log.tran, lm(linkfun(cell.norm) ~ neuron * temp + genotype, data = df)) # F = 21.804, p < 0.001 for gt
# lm5 <- with(log.tran, lm(linkfun(cell.norm) ~ neuron * temp * genotype, data = df)) # no sig interaction
# lm6 <- with(log.tran, lm(linkfun(cell.norm) ~ neuron + temp + genotype, data = df))
# lm7 <- with(log.tran, lm(linkfun(cell.norm) ~ neuron + temp * genotype, data = df))

anova(with(log.tran,lm(linkfun(cell.norm) ~ temp, data = df[df$neuron == "ASI",])),with(log.tran,lm(linkfun(cell.norm) ~ temp+genotype, data = df[df$neuron == "ASI",])))

## Analysis of Variance Table
## 
## Model 1: linkfun(cell.norm) ~ temp
## Model 2: linkfun(cell.norm) ~ temp + genotype
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1     98 4.4439                              
## 2     97 4.2676  1   0.17634 4.0081 0.04807 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# for ASI, F = 4.0071, 1Df p = 0.0481 for genotype
anova(with(log.tran,lm(linkfun(cell.norm) ~ temp, data = df[df$neuron == "ASJ",])),with(log.tran,lm(linkfun(cell.norm) ~ temp+genotype, data = df[df$neuron == "ASJ",])))

## Analysis of Variance Table
## 
## Model 1: linkfun(cell.norm) ~ temp
## Model 2: linkfun(cell.norm) ~ temp + genotype
##   Res.Df    RSS Df Sum of Sq      F        Pr(>F)    
## 1     99 1.6937                                      
## 2     98 1.2197  1   0.47395 38.079 0.00000001533 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# for ASJ, F = 38.082, 1Df, p < 0.001 for genotype
# no significant interaction of temp with genotype for either neuron

lm.ASI <- with(log.tran, lm(linkfun(cell.norm) ~ temp + genotype, data = df[df$neuron == "ASI",]))
lm.ASJ <- with(log.tran, lm(linkfun(cell.norm) ~ temp + genotype, data = df[df$neuron == "ASJ",]))

contrasts.ASI<-summary(lsmeans(lm.ASI, pairwise ~ genotype | temp)$contrasts) %>% dauergut::prange()
contrasts.ASJ<-summary(lsmeans(lm.ASJ, pairwise ~ genotype | temp)$contrasts) %>% dauergut::prange()

stan.ASJ<-with(log.tran, stan_lmer(linkfun(cell.norm) ~ group.id + (1|genoID), data = df[df$neuron == "ASJ",], adapt_delta = 0.99, iter = 4000))

stan.ASI<-with(log.tran, stan_lmer(linkfun(cell.norm) ~ group.id + (1|genoID), data = df[df$neuron == "ASI",], adapt_delta = 0.99, iter = 4000))

mixed.ASI <- stan.ASI %>% getStan_CIs_log_temp(.,group = genotype, base=10)
mixed.ASJ <- stan.ASJ %>% getStan_CIs_log_temp(.,group = genotype, base=10)

#### plot ASI ###
mixed <- mixed.ASI
(p1 <- df %>% dplyr::filter(neuron == "ASI") %>% ggplot(aes(x = genotype, y = cell.norm)) + 
  geom_quasirandom(aes(y=cell.norm, pch = temp),colour = "#339900", cex=1,
                           width = 0.075,size=0.3,
                           method = 'smiley') + 
  list(add.median(width=0.25),add.quartiles()) + facet_grid(.~temp) +
  stat_summary(aes(x=as.numeric(as.factor(genotype)) + 0.3, y=0),
                   fun.data = fun_length, geom = "text", size = 3) +
  add.Bayes.CI() +
  coord_cartesian(ylim=c(0,4095))+
  theme_classic())

#### plot ASJ ###
mixed<-mixed.ASJ
(p2 <- df %>% dplyr::filter(neuron == "ASJ") %>% ggplot(aes(x = genotype, y = cell.norm)) + 
  geom_quasirandom(aes(y=cell.norm, pch = temp),colour = "#339900", cex=1,
                           width = 0.075,size=0.3,
                           method = 'smiley') + 
  list(add.median(width=0.25),add.quartiles()) + facet_grid(.~temp) +
  stat_summary(aes(x=as.numeric(as.factor(genotype)) + 0.3, y=0),
                   fun.data = fun_length, geom = "text", size = 3) +
  add.Bayes.CI() +
  coord_cartesian(ylim=c(0,4095))+
  theme_classic())

__Figure 4 rict-1 are deficient in food suppression

4A

strains<-c("N2", "rict-1(mg360)", "rict-1(ft7)")
#dates<-c("12_5_14", "12_19_14", "1_19_15","3_22_14", "2_25_14", "4_11_14")
dates<-c("12_5_14", "12_19_14","3_22_14", "2_25_14", "4_11_14")
foods <- c("OP50", "HB101")
rict1.food<-read.csv(file.path(pathname, "data", "1B_1C_rict-1_TORC2.csv"), header=TRUE) %>% format_dauer() %>%
  dplyr::filter(day %in% dates) %>% mutate(logit.p = car::logit(pct, adjust=0.01), 
                                           plate.ID = interaction(food, plateID))

# rict1.food %>% ggplot(aes(x=food, y=pct)) + geom_boxplot() + stat_summary(aes(y=pct, group=day, colour = day), fun.y = mean, geom = "line") + facet_grid(.~genotype)

#lm
lm.add <- lm(pct ~ genotype + food, data = rict1.food)
lm.int <- update(lm.add,.~. + genotype*food)
#stan
rict.food.groups <- rict1.food %>% mutate(group.id = interaction(genotype, food))
stan.glmm <- rict.food.groups %>% dauer_stan_groups()

contrasts<-dauergut::dunnett_contrasts(lm.int, ref.index = 1, factor = "genotype", interaction = "food")
mixed<-stan.glmm %>% dauer_StanCIs_groups()

plot.contrasts <- c("",contrasts$interaction$prange[1], "", contrasts$interaction$prange[2], "", contrasts$interaction$prange[3])
plot.contrasts.interaction <- summary(lm.int)$coefficients[,4][5:6] %>% data.frame() %>% mutate(p.value = c(.[[1]],.[[2]]), genotype = strains[2:3]) %>% dauergut::prange()

index<-rep(seq_len(length(strains)), each = 2)

plot.contrasts.H0 <- data.frame(cbind(food = foods, genotype = strains[index])) %>%
  mutate(prange = plot.contrasts)

(p<-ggplot(rict1.food, aes(x=food)) +
  stat_summary(aes(y=pct, group=day), colour = "black", fun.y = mean, geom = "line", alpha = 0.2, linetype = 2) +
  add.median.dauer() +
  add.Bayes.CI() +
  geom_dotplot(aes(y=pct, colour = food, fill=food),binwidth=.015, binaxis="y", position="dodge", stackdir="center", size =.3) +
    #geom_point(aes(y=pct,colour = food), size = 0.7, alpha = 0.75) +
  labs(title = "rict-1 mutants are deficient in food suppression",
           y = "proportion dauer",
           x = "food") +
  facet_grid(.~genotype, switch="both") +
  scale_colour_manual(values = c("black", "#FF9933")) +
    scale_fill_manual(values = c("black", "#FF9933")) +
      scale_y_continuous(breaks=c(0,0.25,0.5,0.75, 1.0)) +
      scale_x_discrete(labels=function(x) sub(" ","\n",x,fixed=TRUE)) +
  geom_text(data = plot.contrasts.H0, aes(x=2, label = prange, y = 1.075, group = NULL), size = 4) +
    stat_summary(aes(x=as.numeric(as.factor(food)) + 0.3, y=-0.05),
                   fun.data = fun_length, geom = "text", size = 3) +
  theme_my_classic + 
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        strip.text.x = element_text(size=10)))

knitr::kable(car::Anova(lm.int), caption="significant interaction of genotype on food")

significant interaction of genotype on food
	Sum Sq	Df	F value	Pr(>F)
genotype	6.3966735	2	233.790975	0.0000000
food	0.0976520	1	7.138133	0.0104016
genotype:food	0.1424647	2	5.206918	0.0091529
Residuals	0.6292950	46	NA	NA

# lsmip(lm.int, genotype ~ food, main="interaction plot")
# 
# plot(rict1.lsm, comparisons = TRUE, alpha = .05, main="confidence intervals") # showing confidence intervals and gap by which they cover zero
# 
# knitr::kable(summary(pairs(rict1.lsm)), caption="contrast matrix for food effect in rict-1 mutants")

_4B #ASJ::daf-28

#days = as.factor(8:10)
strains = c("N2","mg360")
foods = c("OP50","HB101")
days = "1_20_15"

d28 <- read.csv('data/2D_3B_daf-28_GFP_rict-1.csv') %>%
separate(ID, c("ID.A", "ID.B"), sep = ":", extra = "drop") %>%
  subset(food %in% foods & date %in% days & temp == "27" & genotype %in% strains & pheromone == 0 & neuron == "ASJ") %>% 
  mutate(day = as.factor(day), genotype = factor(genotype, levels = strains),
         genoID = interaction(date,genotype, ID.A), food = factor(food, levels = foods), 
         neuron = factor(neuron, levels = "ASJ"), cell.norm = mean)

lm.add <- lm(cell.norm ~ genotype + food, data = d28)
lm.int <- lm(cell.norm ~ genotype * food, data = d28)
anova(lm.add, lm.int) # F = 7.76, Df = 1, p = 0.0077 for genotype interaction

## Analysis of Variance Table
## 
## Model 1: cell.norm ~ genotype + food
## Model 2: cell.norm ~ genotype * food
##   Res.Df      RSS Df Sum of Sq      F   Pr(>F)   
## 1     47 37320917                                
## 2     46 31933165  1   5387751 7.7611 0.007729 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

d28 %<>% mutate(group.id = interaction(genotype, food))

stan.lm <- with(log.tran, stan_glm(formula = linkfun(mean) ~ group.id,
                       data=d28,
                       family = gaussian(),
                       prior = cauchy(),
                       prior_intercept = cauchy(),
                       chains = 3, cores =4, seed = 2000,iter=6000,
                       control = list(adapt_delta=0.99)))

contrasts <- dauergut::dunnett_contrasts(lm.int, ref.index = 1, factor = "genotype", interaction = "food")
mixed <- getStan_CIs_log_groups(stan.lm, group = "food", base = 10) %>% 
  mutate(food = factor(food, levels = foods), genotype = factor(genotype, levels = strains))

#make contrasts with ID for facetting
plot.contrasts <- c("",contrasts$interaction$prange[1],"",contrasts$interaction$prange[2])

index<-rep(seq_len(length(strains)), each = length(foods))
plot.contrasts.H0 <- data.frame(cbind(food = rep(foods, length(strains)), genotype = strains[index])) %>%
  mutate(prange = plot.contrasts)

(p<-ggplot(d28, aes(x=food)) +
  #stat_summary(aes(y=mean, group=day), colour = "black", fun.y = mean, geom = "line", alpha = 0.2, linetype = 2) +
  list(add.median(),add.quartiles()) +
       geom_quasirandom(aes(y=cell.norm, ),colour = "#339900", cex=1,
                           width = 0.075,size=0.3,
                           method = 'smiley') +
  #geom_point(aes(y=cell.norm,colour = food), size = 0.7, alpha = 0.75) +
  labs(title = "rict-1 mutants are deficient in food suppression",
           y = "mean GFP expression",
           x = "food") +
    geom_text(data = plot.contrasts.H0, aes(x=food, label = prange, y = 4200, group = NULL), size = 4) + 
  facet_grid(.~genotype, switch="both") +
  scale_colour_manual(values = c("black", "#FF9933")) +
  theme_my_classic +
    stat_summary(aes(x=as.numeric(as.factor(food)) + 0.3, y=100),
                   fun.data = fun_length, geom = "text", size = 3))

Fig 5 intestinal mTORC2 alters foraging behaviors

5A

intestinal rict-1 expression rescues behavior.

5B

roam<-read.csv(file.path(pathname, "data/5B_5E_gridentry.csv")) %>% 
  mutate(strain = interaction(genotype, line, drop=TRUE),
         strainDate = interaction(strain, date, drop=TRUE),
         plateID = interaction(strainDate,plate, drop=TRUE),
         total.boxes = 186)
#reorder to set N2 = level1
# newlevels<-subset(roam, genotype =! "N2")
# newlevels<-subset(levels(newdata$genotype))
# 
# newlevels<-newlevels[newlevels$genotype, drop=TRUE]
# newlevels$genotype <- factor(newlevels$genotype, drop=TRUE)
# roam$genotype <- factor(roam$genotype, levels = c("N2", levels(newlevels$genotype))

library(lme4)
#simple ANOVA:
lms<-list(lm = lm(n_entries ~ genotype, data = roam),
          lm.log = lm(log(n_entries) ~ genotype, data = roam),
          lm.nboxes = lm(n_boxes ~ genotype, data = roam))

#glmms
glmms<-list(glmm.date = glmer(data=roam, n_entries ~ genotype + (1|date), family="poisson"),
            glmm.strainDate = glmer(data=roam, n_entries ~ genotype + (1|strainDate), family="poisson"),
            glmm.obs = glmer(data=roam, n_entries ~ genotype + (1|date) + (1|plateID), family="poisson"),
            glmm.obsOnly = glmer(data=roam, n_entries ~ genotype + (1|plateID), family="poisson"),
            glmm.fixDate = glmer(data=roam, n_entries ~ date + genotype + (1|plateID), family="poisson"),
            glmm.fixDate.ranStrain = glmer(data=roam, n_entries ~ date + genotype + (1|plateID) + (1|strain), family="poisson"),
            glmm.nest = glmer(data=roam, n_entries ~ genotype + (1|date/plateID), family="poisson"), #nested individual random effect
            glmm.binom = glmer(data=roam,
                               formula = cbind(roam$n_boxes, roam$total.boxes - roam$n_boxes) ~ genotype + (1|date/plateID),
                               family="binomial"))

#non mixed glms
glms <- list(glm = glm(data=roam, n_entries ~ date + genotype, family="poisson"),
             glm.nb = glm.nb(data=roam, n_entries ~ date + genotype, init.theta = 1, link="log"),
             glm.binom = glm(data=roam, cbind(roam$n_boxes, roam$total.boxes - roam$n_boxes) ~ date + genotype, family="binomial"))

mods<-c(lms, glms, glmms)


###checking distributional assumptions:
# library(DHARMa)

# sim.resids<-lapply(mods[1:length(mods)], function(x) {
#   simulationOutput <- simulateResiduals(fittedModel = x, n=250, refit = T)
#   #par(mfrow = c(2,2))
#   #plotResiduals(roam$genotype, simulationOutput$scaledResiduals)
#   #plotResiduals(roam$date, simulationOutput$scaledResiduals)
#   #plotResiduals(roam$strain, simulationOutput$scaledResiduals)
#   return(simulationOutput)
#   })
# 
# lapply(sim.resids, function(x) {
#   plotSimulatedResiduals(simulationOutput = x)
# })


# using DHARMa package check distribution of residuals by predictors (takes a while using refit = T):
# (mods.dispersion <-
#   lapply(glmms, function(x) {
#   simulationOutput <- simulateResiduals(fittedModel = x, n=250, refit = T)
#   testOverdispersion(simulationOutput = simulationOutput)
# }))

# nested glmm looks rougly the best for both count date (n_entries) and binary data (n_boxes/total) 
# although there is still some small effect of 
# predicted value on residual quantiles

strains <- c("N2", "ft7", "ft7;ex[ges-1p]", "ft7;ex[rict-1p]")
days <- c("10_1_16", "10_5_16", "10_9_16")
## remove outliers:

roam %<>% dauergut::flag_outliers(df = ., lin.mod = mods$lm, threshold = 4)

roam.resc <- roam %>% subset(genotype %in% strains & date %in% days) %>% 
  mutate(genotype = factor(genotype, levels = strains))

glmm.nest <- update(mods$glmm.nest, data = roam.resc[roam.resc$outlier.status == FALSE,] )
stan.glmm <- stan_glmer(data = roam.resc[roam.resc$outlier.status == FALSE,], 
                        formula = n_entries ~ genotype + (1|date/plateID), family = "poisson",
                        iter = 4000)
  
contrasts<-dauergut::tukey_contrasts(glmm.nest, "genotype")
mixed<-stan.glmm %>% roam_getStan_CIs()
plot.contrasts<-c("",contrasts$prange[1], "","")
plot.contrasts.2<-c("","",contrasts$prange[4:5])



(p<-dauergut::plot_CIs(roam.resc, title='mTORC2 components acts in the intestine to regulate foraging', 
                       plot.contrasts=plot.contrasts, plot.contrasts.2 = plot.contrasts.2,ypos = 800, offset = 0, type = "grid"))

# #print this at the end
# for (i in sim.resids) {
#   print(summary(i$fittedModel))
#   plotSimulatedResiduals(simulationOutput = i)
# }

S2A

strains<-c("N2","ft7")
d7GFP<-read.csv("data/S2A_daf7GFP_adults.csv") %>% filter(mean!=4095 & genotype %in% strains) %>% mutate(genotype = factor(genotype, levels=strains), ID = as.character(ID)) %>%
  separate(ID, c("ID.A", "ID.B","ID.C"), sep = ":", extra = "drop") %>% 
  mutate(genoID = paste(genotype, ID.A, animal, sep = ":")) #genoID is animal, 2 cells per animal measured.

df <- d7GFP %>% group_by(genotype, genoID) %>% summarise(cell.norm = mean(mean)) # take mean of each worm
library(lsmeans)
# log.tran<-make.tran(type="genlog", param = c(0.5,10))
log.tran<-make.tran(type="genlog", param = c(0,10))
### t-test with log-transformed data
d7.t <-with(log.tran,(t.test(data=df,linkfun(cell.norm)~genotype)))

N2<-df[df$genotype=="N2","cell.norm"] %>% as.matrix
ft7<-df[df$genotype=="ft7","cell.norm"] %>% as.matrix
BEST.t <- with(log.tran, BEST::BESTmcmc(linkfun(N2),linkfun(ft7))) # variance is equal effects sizes are same as stan_glm

## 
## Processing function input....... 
## 
## Done. 
##  
## Beginning parallel processing using 3 cores. Console output will be suppressed.
## 
## Parallel processing completed.
## 
## MCMC took 0.167 minutes.

stanlm <- with(log.tran, stan_glm(linkfun(cell.norm) ~ genotype, family=gaussian(), data = df, prior = cauchy(), prior_intercept = cauchy()))

contrasts <-data.frame(p.value = d7.t$p.value) %>% dauergut::prange()
plot.contrasts <- c("",contrasts$prange)
mixed <- stanlm %>% getStan_CIs_log(base = 10)

(p<-plot_CIs(df, plot.contrasts = plot.contrasts, type = "GFP", title = "", ypos = 4700, offset = 0))

5C

#### see speed_curvature function for data analysis

rict-1 mutation increases roaming behavior

5D

devtools::install_github("rasmusab/bayesian_first_aid")
strains <- c("N2","ft7")
roam_beh <- read.csv(file.path(pathname, "data/5D_roam_dwell_all.csv")) %>% 
  mutate(genotype = factor(genotype, levels = strains), pct = pct.roam)


roam.t <- t.test(roam_beh[roam_beh$genotyp == "N2",]$pct, roam_beh[roam_beh$genotyp == "ft7",]$pct)
BEST.t <- BEST::BESTmcmc(roam_beh[roam_beh$genotyp == "N2",]$pct, roam_beh[roam_beh$genotyp == "ft7",]$pct)

## 
## Processing function input....... 
## 
## Done. 
##  
## Beginning parallel processing using 3 cores. Console output will be suppressed.
## 
## Parallel processing completed.
## 
## MCMC took 0.049 minutes.

summary <- lapply(strains, function(x) {
  summary <- roam_beh %>% dplyr::filter(genotype == x) %>% 
  summarize(roam = sum(roam), n = (sum(roam,dwell)))
}) %>% unlist()

roamers <- summary[c(1,3)]
total <- summary[c(2,4)]

prop.test(roamers, total)

## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  roamers out of total
## X-squared = 2156.2, df = 1, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.2376133 -0.2181689
## sample estimates:
##    prop 1    prop 2 
## 0.1177278 0.3456189

fit <- BayesianFirstAid::bayes.prop.test(roamers, total)

glmm <- lme4::glmer(data = roam_beh, cbind(roam, dwell) ~ genotype + (1|date:plate), family = binomial)

stan.glmm <- stan_glmer(data = roam_beh, 
                        formula = cbind(roam, dwell) ~ genotype + (1|date:plate), 
                        family = binomial,
                        iter = 4000)

mixed <- dauer_getStan_CIs(stan.glmm)
contrasts <- data.frame(p.value = roam.t$p.value) %>% dauergut::prange()
plot.contrasts <- c("",contrasts$prange)


(p<-dauergut::plot_CIs(roam_beh, plot.contrasts = plot.contrasts, type = "dauer", title = "", ypos = 1.05, offset = 0) + labs(subtitle = "", y = "proportion roaming"))

5E

pdfr-1 is epistatic to rict-1 in foraging behavior.

  roam<-read.csv(file.path(pathname, "data/5B_5E_gridentry.csv")) %>% 
  mutate(strain = interaction(genotype, line, drop=TRUE),
         strainDate = interaction(strain, date, drop=TRUE),
         plateID = interaction(strainDate,plate, drop=TRUE),
         total.boxes = 186)

strains <- c("N2", "ft7", "pdf-1", "pdfr-1", "pdfr-1;ft7")
days <- c("12_20_16", "1_5_17", "1_8_17")
## remove outliers (based on analysis, below):
lm <- lm(n_entries ~ genotype, data = roam)
roam %<>% dauergut::flag_outliers(df = ., lin.mod = lm, threshold = 4)

roam.pdfr1 <- roam %>% subset(genotype %in% strains & date %in% days) %>%
  mutate(genotype = factor(genotype, levels = strains))

glmm.nest <-glmer(data = roam.pdfr1[roam.pdfr1$outlier.status == FALSE,], formula = n_entries ~ genotype + (1|date/plateID), family = "poisson")
stan.glmm <- stan_glmer(data = roam.pdfr1[roam.pdfr1$outlier.status == FALSE,], 
                        formula = n_entries ~ genotype + (1|date/plateID), family = "poisson",
                        iter = 4000, adapt_delta = 0.99)

contrasts<-dauergut::tukey_contrasts(glmm.nest, "genotype")
mixed<-stan.glmm %>% roam_getStan_CIs()
plot.contrasts<-c("",contrasts$prange[1:4])
plot.contrasts.2<-c("","","","",contrasts$prange[10])

(p<-dauergut::plot_CIs(roam.pdfr1, title='rict-1 effects on roaming behavior depends on pdf-1/pdfr-1 signaling', plot.contrasts=plot.contrasts, plot.contrasts.2 = plot.contrasts.2,ypos = 800, offset = 50, type = "grid"))

5F

pdf-1::venus is increased in both gut and coelemocytes in rict-1 mg360

Fig 6 Natural variation in thermally-driven dauer formation

6A

#knitr::opts_chunk$set(cache=TRUE)
filepath<-file.path(pathname, "data/6A_wildHTdata_allbinary.csv")
strains <- c("N2","CB4856","JU561", "MY14","DL238","CB4852","QX1211","JU322","JU775","JU1400","ED3072","AB1","AB3","JU362","ED3040", "JU345","PX178","CB3198")

##"TR403","KR314",CB4507,DR1350 left off due to low n
HT.wildall<-read.csv(filepath, header=TRUE) %>% format_dauer()

lm <- HT.wildall %>% dauer_ANOVA()
contrasts<-dunnett_contrasts(lm, ref.index = 1, "genotype")

## [1] "interaction term not indicated"

stan.glmm <- HT.wildall %>% dauer_stan_glmm()
mixed<-stan.glmm %>% dauer_getStan_CIs()

plot.contrasts<-c("",contrasts$prange[1:17])

(p<-plot_CIs(HT.wildall, title='C. elegans wild isolates vary in temp dependent dauer formation', plot.contrasts, ypos = 1.075, type = "dauer"))

#plot

Wild isolates were tested for sensitivity to high temperature-induced dauer formation at 27ºC. Most strains were tested in duplicate over 3 biologically independent assays.

S3A-B

filepath<-file.path(pathname, "data/S4B_C3_sumstats.csv")
C327<-read.csv(filepath, header=TRUE)

(p<-dauergut::plot_regression_se(C327, "mean.C3", "mean.HT","SEM.C3", "SEM.HT"))

6B

strains <- c("N2","CB4856", "JU561", "MY14", "JU322", "JU362","ED3040", "JU345", "PX178", "CB3198")

wild.roam<-read.csv(file.path(pathname, "data/6B_7F_roam_wild_isolates.csv")) %>%
  dplyr::filter(genotype %in% strains) %>%
  mutate(genotype = factor(genotype, levels = strains),
         strainDate = interaction(genotype, date, drop=TRUE),
         plateID = interaction(strainDate,plate, drop=TRUE),
         total.boxes = 186)

lm <- lm(n_entries ~ genotype, data = wild.roam)

wild.roam %<>% dauergut::flag_outliers(df = ., lin.mod = lm, threshold = 4)

stan.glmm <- stan_glmer(data = wild.roam[wild.roam$outlier.status == FALSE,], 
                        formula = n_entries ~ genotype + (1|date/plateID), family = "poisson",
                        iter = 4000, adapt_delta = 0.99)

lm <- update(lm, data = wild.roam[wild.roam$outlier.status == FALSE,])

contrasts<-dauergut::tukey_contrasts(lm, "genotype")
mixed<-stan.glmm %>% roam_getStan_CIs()
plot.contrasts<-c("",contrasts$prange[1:9])

(p<-plot_CIs(wild.roam, title='Foraging behavior is not correlated with Hid phenotypes', plot.contrasts=plot.contrasts, ypos = 800, offset = 50, type = "grid"))

Fig 7 Identification of two QTL for 27º dauer formation

7A

filepath<-file.path(pathname,'data/7A_CSSdata_allbinary.csv')
strains<- c("N2","CB4856","CSSI","CSSII","CSSIII","CSSIV","CSSV","CSSX")
CSS<-read.csv(filepath, header=TRUE) %>% format_dauer()

lm <- CSS %>% dauer_ANOVA()
glmm.locus <- glmer(CSS,formula = cbind(dauer,non) ~ gtChrI + gtChrII + gtChrIII + gtChrIV + gtChrV + gtChrX + (1|plateID) + (1|strainDate), family = binomial(link = logit),control=glmerControl(optimizer="bobyqa"))
contrasts<-dauergut::dunnett_contrasts(lm, ref.index = 1, "genotype")

## [1] "interaction term not indicated"

plot.contrasts<-c("",contrasts$prange[1:7])

stan.glmm <- CSS %>% dauer_stan_glmm()
mixed<- stan.glmm %>% dauer_getStan_CIs() #CIs

(p<-dauergut::plot_CIs(CSS, title='CB4856 Chromosome II contains 1 or more dauer QTLs', plot.contrasts, ypos = 1.075, type = "dauer"))

knitr::kable(drop1(glmm.locus, test = "Chisq"))

	Df	AIC	LRT	Pr(Chi)
	NA	263.5956	NA	NA
gtChrI	1	261.7947	0.1991198	0.6554323
gtChrII	1	277.8591	16.2635395	0.0000551
gtChrIII	1	262.3102	0.7146253	0.3979126
gtChrIV	1	264.3476	2.7519800	0.0971341
gtChrV	1	262.3788	0.7831660	0.3761746
gtChrX	1	265.0531	3.4575485	0.0629636

6A. Chromosome substitution strains demonstrate a chromosome II QTL for 27º dauer formation. Chromosome-substitution strains (each containing a single CB4856 chromosome introgressed into N2) were assayed for dauer formation at 27º. Proportion of dauers in each strain is shown. Box and scatter plot show raw data. Bayesian 90% (grey) and 50% (red) credible intervals are shown on the right (see methods). All p-values reflect ANOVA with Dunnett post-hoc comparision to N2. All p-values reflect Dunnett post-hoc comparision to N2. n=6 experiments over 3 independent days. Likelihood ratio test for genotype-by-locus shown below.

7C

filepath<-file.path(pathname, "data/7C_ewIRdata_allbinary.csv")

strains=c("N2","CB4856","ewIR19","ewIR20",
          "ewIR21","ewIR22","ewIR23","ewIR24",
          "ewIR25","ewIR26","ewIR27")

ewIR<-read.csv(filepath, header=TRUE) %>% format_dauer()


ewIR.nocont<-subset(ewIR, !(genotype %in% c("N2","CB4856")))
mapping<-read.csv(file.path(pathname, 'data/7C_chrII.markers.concise.csv'))
mapping.nocont<-subset(mapping, !(genotype %in% c("N2","CB4856")))
mapping.nocont$other.loci<-NULL
ewIR<-merge(ewIR, mapping, by="genotype", all=TRUE)
ewIR.nocont<-merge(ewIR.nocont, mapping.nocont, by="genotype", all=TRUE)

#glmm.nest = glmer(data = ewIR, cbind(dauer,non) ~ genotype + (1|day/plateID), family="binomial",control=glmerControl(optimizer="bobyqa"))
lm <- ewIR %>% dauer_ANOVA()
contrasts<-dauergut::tukey_contrasts(lm, "genotype")
stan.glmm <- ewIR %>% dauer_stan_glmm()

mixed<- stan.glmm %>% dauer_getStan_CIs()
plot.contrasts<-c("",contrasts$prange[1:10])
# get contrasts compared to ewIR25
#plot.contrasts.2<-c("","",contrasts$prange[c(25,32,38,43,47,50)], "", contrasts$prange[53:54]) 
plot.contrasts.2 <- c("","","","","","","","","",contrasts$prange[53:54])


(p<-dauergut::plot_CIs(ewIR, title='CB4856 27º dauer formation QTL lies on ChrII', plot.contrasts, plot.contrasts.2 = plot.contrasts.2,ypos = 1.075, offset = 0.075, type = "dauer"))

knitr::include_graphics('/Users/mikeod/Dropbox/Temp_shareDocs/Dauer/rict-1_paper/Temperature_version/fig2/ewIR-map.png')

7B. Chromosome near isogenic lines (NILs) were tested for 27º dauer formation. Each strain contains a segment of CB4856 DNA (blue, right) introgressed into an N2 background. EwIR25 contains a QTL which reduces dauer formation. ewIR27 contains an additional QTL that increases dauer formation. Diagram shows putative location of these QTL for decreased (blue) or increased (red) dauer formation. Box and scatter plot show raw data. Blue dot and bars show marginal means and 95% confidence interval, respectively, generalized linear mixed effects model using strain as a predictor (see methods). Black p-values reflect Tukey post-hoc comparision to N2, red p-values indicate comparison to ewIR25. n≥6 assays on ≥3 separate days.

7E

file<-'7E_MONIL_allbinary.csv'
filepath<-file.path(pathname, "data", file)
strains<-c("N2","CB4856","NIL10","NIL59",
                               "NIL76","NIL78")

MONIL<-read.csv(filepath, header=TRUE) %>% format_dauer()
  # mutate(genotype = factor(genotype, levels = strains),
  #        pct = as.numeric(paste(dauer/(dauer+non))),
  #        strainDate = interaction(genotype, day),
  #        plateID = interaction(strainDate, plate))

file2<-'7E_MONIL_mapping.csv'
filepath2<-file.path(pathname, "data", file2)
mapping<-read.csv(filepath2)
controls <- c("N2", "CB")

MONIL <- merge(MONIL, mapping, by="genotype", all=TRUE) %>% mutate(binII12680000 = factor(binII12680000, levels = controls), 
                 binII12790000 = factor(binII12790000, levels = controls),
                 binII12900000 = factor(binII12900000, levels = controls),
                 binII13050000 = factor(binII13050000, levels = controls),
                 other = factor(other, levels = controls))

glmm.nest = glmer(data = MONIL, cbind(dauer,non) ~ genotype + (1|day/plateID), family="binomial", control=glmerControl(optimizer="bobyqa"))

lm <- MONIL %>% dauer_ANOVA()
stanlmer <- MONIL %>% dauer_stan_glmm()
#per-locus glmm - loci are fixed effects
loci<-paste(names(mapping[-1]), collapse=" + ")
rand<-"(1|day/plateID)"
f<-paste("cbind(dauer,non) ~",loci," + ",rand)
f.rev <- paste("cbind(dauer,non) ~",rev(loci)," + ",rand)
locus.glmm<-glmer(f, data=MONIL, family=binomial, control=glmerControl(optimizer="bobyqa"))
locus.rev.glmm <-glmer(f.rev, data=MONIL, family=binomial, control=glmerControl(optimizer="bobyqa"))

# #bootstrapping
# # Not run:
# nsim=100 # takes ~12 seconds per iteration
# PBmod<-drop1(locus.glmer,test="user",sumFun=PBSumFun)
# PBmod
# # end (**Not run**)

contrasts<-dauergut::tukey_contrasts(glmm.nest, "genotype")
mixed <- dauer_getStan_CIs(stanlmer)
plot.contrasts<-c("",contrasts$prange[1:5])
(p<-plot_CIs(MONIL, title='A ~100kb interval defines QTL1', plot.contrasts, ypos = 0.8, type = "dauer"))

7F

strains<- c("N2","CB4856","QTL1","QTL2")
days <- c("3_10_17", "3_17_17", "4_6_17")
QTL.roam<-read.csv(file.path(pathname, "data/6B_7F_roam_wild_isolates.csv")) %>%
  dplyr::filter(genotype %in% strains & date %in% days) %>%
  mutate(genotype = factor(genotype, levels = strains),
         strainDate = interaction(genotype, date, drop=TRUE),
         plateID = interaction(strainDate,plate, drop=TRUE),
         total.boxes = 186)


lm <- lm(n_entries ~ genotype, data = QTL.roam)

QTL.roam %<>% dauergut::flag_outliers(df = ., lin.mod = lm, threshold = 4)

glmm.nest = glmer(data=QTL.roam, n_entries ~ genotype + (1|date/plateID), family="poisson")

stan.glmm <- stan_glmer(data = QTL.roam[QTL.roam$outlier.status == FALSE,], 
                        formula = n_entries ~ genotype + (1|date/plateID), family = "poisson",
                        iter = 4000, adapt_delta = 0.99)

lm <- update(lm, data = QTL.roam[QTL.roam$outlier.status == FALSE,])

contrasts<-dauergut::tukey_contrasts(lm, "genotype")
mixed<-stan.glmm %>% roam_getStan_CIs()
plot.contrasts<-c("",contrasts$prange[1:3])
plot.contrasts.2 <- c("","",contrasts$prange[4:5])

(p<-dauergut::plot_CIs(QTL.roam, title='Foraging behavior is not correlated with Hid phenotypes', plot.contrasts=plot.contrasts, plot.contrasts.2 = plot.contrasts.2, ypos = 800, offset = 50, type = "grid"))

contrasts<-dauergut::tukey_contrasts(glmm.nest, factor="genotype")

plot.contrasts<-c("",contrasts$prange[1:3])
plot.contrasts.2<-c("","", contrasts$prange[4:5])

(P<-plot_CIs(QTL.roam, title='Chromosome II may encode a foraging QTL', plot.contrasts=plot.contrasts, plot.contrasts.2 = plot.contrasts.2, ypos = 800, offset = 50, type = "grid"))

__Supplemental fig 4 - analysis of QTL1 breakpoints

S4C

strains<-c("N2", "ewIR25", "NK", "NK_DG")
foods <- "OP50"

NKDG<-read.csv(file.path(pathname, "data", "S4C_NKDG_allbinary.csv"), header=TRUE) %>% format_dauer_pd() #including partial/pd as non-dauer due to excess zeros. 

##glmm
#glmm <- TORC2 %>% dauer_pd_glmm() #including partial/pd as non-dauer due to excess zeros. 
#lm
lm <- NKDG %>% dauer_ANOVA()
#stan
stan.glmm <- NKDG %>% dauer_stan_glmm()

contrasts<-dauergut::dunnett_contrasts(lm, ref.index = 1, "genotype")

## [1] "interaction term not indicated"

mixed<-stan.glmm %>% dauer_getStan_CIs()
plot.contrasts<-c("",contrasts$prange[1:3])

(p<-dauergut::plot_CIs(NKDG, title='mTORC2 prevents high temperature dauer formation', plot.contrasts, ypos = 1.075, type = "dauer"))

Dauer intestine rict-1 Paper

Mike O’Donnell